Package “"proteoQ"”

Qiang Zhang

2019-05-21

Introduction to proteoQ

Chemical labeling using tandem mass tag (TMT) has been commonly applied in mass spectrometry (MS)-based quantification of proteins and peptides. The proteoQ tool is designed to aid automated and reproducible analysis of proteomics data. It interacts with an Excel spread sheet for dynamic sample selections, aesthetics controls and statistical modelings. The arrangement allows users to put data manipulation behind the scene and apply metadata to quickly address interesting biological questions using various informatic tools. In addition, the entire workflow is documented and can be conveniently reproduced upon revisiting.

The tool currently processes the peptide spectrum matches (PSM) tables from Mascot searches for 6-, 10- or 11-plex TMT experiments. Peptide and protein results are then produced with users’ selection of parameters in data filtration, alignment and normalization. The package further offers a suite of tools and functionalities in statistics, informatics and data visualization by creating ‘wrappers’ around published R functions.

Installation

To install this package, start R (version “3.6”) and enter:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("Biobase", "GSVA", "Mfuzz", "gage", "limma"))

if (!requireNamespace("devtools", quietly = TRUE))
    install.packages("devtools")
devtools::install_github("qzhang503/proteoQ@master")

Application

In this section I illustrate the following applications of proteoQ:

Set up the experiments

PSM table(s) in a csv format will be exported by the users from the Mascot search engine. I typically set the option of Include sub-set protein hits to 0 with my opinionated choice in satisfying the principle of parsimony. The options of Header and Peptide quantitation should be checked to include the search parameters and quantitative values. The filename(s) of the export(s) will be taken as is.1

The same peptide sequence under different PSM files can be assigned to different protein IDs when inferring proteins from peptides using algorithms such as greedy set cover. To avoid such ambiguity in protein inference, I typically enable the option of Merge MS/MS files into single search in Mascot Daemon. If the option is disabled, peptide sequences that have been assigned to multiple protein IDs will be removed for now when constructing peptide reports.

The pacakge reads an Excel template containing the metadata of multiplex experiment numbers, including TMT channels, LC/MS injection indices, sample IDs, corresponding RAW data file names and addditional fields from the users. The default file name for the experimental summary is expt_smry.xlsx. If samples were fractionated off-line prior to LC/MS, a second Excel template will also be filled out to link multiple RAW file names that are associated to the same sample IDs. The default file name for the fractionation summary is frac_smry.xlsx. The function, extract_raws, can be used to summarise .raw file names under a file folder:

Note that the above files should be stored immediately under the the file folder specified by dat_dir. Examples of PSM outputs, expt_smry and frac_smry can be found as the follows:

and the description of the column keys in the Excel files can be found from the help document:

As a final step of the setup, we will load the experimental summary and some precomputed results:

Summarize PSMs to peptides and proteins

Process PSMs — In this section, I demonstrate the summarisation of PSM data to peptides and proteins. The data set I use in this section corresponds to the proteomics data from Mertins et al.(2018). In the study, two different breast cancer subtypes, WHIM2 and WHIM16, from patient-derived xenograft models were assessed by three independent laborotories. Under each location, lysates from WHIM2 and WHIM16 were each split and labeled with 10-plex TMT at equal sample sizes and repeated on a different day. We start by processing PSM data from Mascot outputs:

PSM outliers will be assessed at a basis of per peptide and per sample at rm_outliers = TRUE, which can be a slow process for large data sets. To circumvent repeated efforts in the assessment of PSM outliers, we may set rm_outliers = FALSE and plot_violins = TRUE when first executing normPSM(). We then visually inspect the violin plots of reporter-ion intensity. Empirically, PSMs with reporter-ion intensity less than 1,000 are trimmed and samples with median intensity that is 2/3 or less to the average of majority samples are removed from further analysis.2

Summarize PSMs to peptides — We next summarise PSM to peptides.

At id = pep_seq_mod, peptide sequences that are different in variable modificaitons will be treated as different species. The log2FC of peptide data will be aligned by median centering across samples by default. If method_align = MGKernel is chosen, log2FC will be aligned under the assumption of multiple Gaussian kernels.3 The parameter n_comp defines the number of Gaussian kernels and seed set a seed for reproducible fittings. The parameters range_log2r and range_int define the range of log2FC and the range of reporter-ion intensity, respectively, for use in the scaling of standard deviation across samples.

Let’s compare the log2FC profiles with and without scaling normalization:4

There are 60 panels of of histograms in each plot, which may not be easy to explore as a whole. In stead, we will break the plots down by their data origins. We begin with modifying the expt_smry.xlsx file by adding the columns BI, JHU and PNNL. Each of the new columns includes sample entries that are tied to their laboratory origins.

Select subsets

We now are ready to plot histograms for each subset of data.5 In the tutorial, we only display the plots using the BI subset:

*NB*: We interactively told `pepHist()` that we are interested in sample entries under the newly created `BI` column. We also supply a file name, assuming that we want to keep the earlierly generated plots with default file names of `Peptide_Histogram_N.png` and `Peptide_Histogram_Z.png`. 
**Figure 1.** Histograms of peptide log2FC. Left: `scale_log2r = FALSE`; right, `scale_log2r = TRUE`**Figure 1.** Histograms of peptide log2FC. Left: `scale_log2r = FALSE`; right, `scale_log2r = TRUE`

Figure 1. Histograms of peptide log2FC. Left: scale_log2r = FALSE; right, scale_log2r = TRUE

As expected, the widths of log2FC profiles become more consistent after the scaling normalization. However, such adjustment may cause artifacts when the standard deviaiton across samples are genuinely different. I typically test scale_log2r at both TRUE and FALSE, then make a choice in data scaling together with my a priori knowledge of the characteristics of samples.6 Alignment of log2FC against housekeeping or normalizer protein(s) is also available. This seems suitable when the quantities of proteins of interest are different across samples where the assumption of constitutive expression for the vast majority of proteins may not hold.

Summarize peptides to proteins — We then summarise peptides to proteins using a two-component Gaussian kernel.

Similar to the peptide summary, we inspect the alignment and the scaling of ratio profiles, and re-normalize the data if needed.7

MDS and PCA plots

In this section, we visualize MDS, PCA and Euclidean distance against the peptide data at scale_log2r = TRUE. We start with metric MDS for peptide data:

**Figure 2A.** MDS of peptide log2FC at `scale_log2r = TRUE`

Figure 2A. MDS of peptide log2FC at scale_log2r = TRUE

It is clear that the WHIM2 and WHIM16 samples are well separated by the Euclidean distance of log2FC (Figure 2A). We next take the JHU data subset as an example to explore batch effects in the proteomic sample handling:

**Figure 2B-2C.** MDS of peptide log2FC for the `JHU` subset. Left: original aesthetics; right, modefied aesthetics**Figure 2B-2C.** MDS of peptide log2FC for the `JHU` subset. Left: original aesthetics; right, modefied aesthetics

Figure 2B-2C. MDS of peptide log2FC for the JHU subset. Left: original aesthetics; right, modefied aesthetics

We immediately spot that all samples are coded with the same color (Figure 2B). This is not a surprise as the values under column expt_smry.xlsx::Color are exclusively JHU for the Select_JHU subset. For similar reasons, the two different batches of TMT1 and TMT2 are distinguished by transparency, which is governed by column expt_smry.xlsx::Alpha. We may wish to modify the aesthetics using different keys: e.g., color coding by WHIMs and size coding by batches, without the recourse of writing new R scripts. One solution is to link the attributes and sample IDs by creating additional columns in expt_smry.xlsx. In this example, we have had coincidentally prepared the column Shape and Alpha to code WHIMs and batches, respectively. Therefore, we can recycle them directly to make a new plot (Figure 2C):

The prnMDS performs MDS for protein data. For PCA analysis, the corresponding functions are pepPCA and prnPCA for peptide and protein data, respectively.

While MDS approximates Euclidean distances at a low dimensional space. Sometime it may be useful to have an accurate view of the distance matrix. Functions pepEucDist and prnEucDist plot the heat maps of Euclidean distance matrix for peptides and proteins, respectively. They are wrappers of (pheatmap) and inherit many parameters therein. Supposed that we are interested in visualizing the distance matrix for the JHU subset:

Parameter annot_cols defines the tracks to be displayed on the top of distrance-matrix plots. In this example, we have choosen expt_smry.xlsx::Shape and expt_smry.xlsx::Alpha, which encodes the WHIM subtypes and the batch numbers, respectively. Parameter annot_colnames allows us to rename the tracks from Shape and Alpha to WHIM and Batch, respectively, for better intuition.

**Figure 2D.** EucDist of peptide log2FC at `scale_log2r = TRUE`

Figure 2D. EucDist of peptide log2FC at scale_log2r = TRUE

Correlation plots

In this section, we visualize the batch effects through correlation plots. The proteoQ tool currently limits itself to a maximum of 44 samples for a correlation plot. In the demo, we will perform correlation analysis against the PNNL data subset. By default, samples will be arranged diagnoally from upper left to bottom right by the row order of expt_smry.xlsx::Sample_ID within a subset. We have learned from the earlier MDS analysis that the batch effects are smaller than the differences between W2 and W16. We may wish to put the TMT1 and TMT2 groups adjacient to each other for visualization of more nuance batch effects, followed by the correlational comparison of WHIM subtypes. We can achieve this by supervising sample IDs at a customized order. In the expt_smry.xlsx, I have prepared an Order column where samples within the JHU subset were arranged in the order of W2.TMT1, W2.TMT2, W16.TMT1 and W16.TMT2. Now we tell the program to look for the Order column for sample ordering:

**Figure 3A-3B.** Correlation of log2FC for the `PNNL` subset. Left: peptide; right, protein**Figure 3A-3B.** Correlation of log2FC for the `PNNL` subset. Left: peptide; right, protein

Figure 3A-3B. Correlation of log2FC for the PNNL subset. Left: peptide; right, protein

The documentation from this point on remains under construction…

The following performs of heat map visualization against protein data:

**Figure 4.** Heat map visualization of protein log2FC at `scale_log2r = TRUE`

Figure 4. Heat map visualization of protein log2FC at scale_log2r = TRUE

The following performs the imputation of peptide and protein data:

The following performs the trend analysis against protein expressions:

The following performs the NMF analysis against protein data:

The following performs the significance analysis of peptide and protein data:

**Figure 5A-5B.** Volcano plot visualization of log2FC. Left: contrasts defined in `contrasts_1`; right, contrasts defined in `contrasts_2`**Figure 5A-5B.** Volcano plot visualization of log2FC. Left: contrasts defined in `contrasts_1`; right, contrasts defined in `contrasts_2`

Figure 5A-5B. Volcano plot visualization of log2FC. Left: contrasts defined in contrasts_1; right, contrasts defined in contrasts_2

The following performs GSVA:

The following maps gene sets under the environment of volcano plot visualization:

Philipp, Martins. 2018. “Reproducible Workflow for Multiplexed Deep-Scale Proteome and Phosphoproteome Analysis of Tumor Tissues by Liquid Chromatography-Mass Spectrometry.” Nature Protocols 13 (7): 1632–61. https://doi.org/10.1038/s41596-018-0006-9.


  1. The default file names begin with letter F, followed by six digits and ends with .csv in file name extension.

  2. The sample removal and PSM re-processing can be achieved by deleting the corresponding entries under the column Sample_ID in expt_smry.xlsx, followed by the re-load of the experiment, load_expts(), and the re-execution of normPSM() with desired parameters.

  3. Density kernel estimates can occasionally capture spikes in the profiles of log2FC for data alignment. Users will need to inspect the alignment of ratio histograms and may optimize the data normalization with different combinations of tuning parameters before proceeding to the next steps.

  4. normPep() will report log2FC results both before and after the scaling of standard deviations.

  5. system files will be automatically updated from the modified expt_smry.xlsx

  6. The default is scale_log2r = TRUE throughout the package. When calling functions involved parameter scale_log2r, users will specify explicitly scale_log2r = FALSE to overwrite the default. Although the package provides the facility to look for a global setting of scale_log2, I don’t recommend using it.

  7. Prameter fasta is solely used for the calculation of protein percent coverage. Precomputed data will be used if no fasta database is provided.